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We propose a theoretical/computational protocol based on the use of the Ground State (GS) Path Integral (PI) 
Quantum Monte Carlo (QMC) for the calculation of the kinetic and Coulomb energy density for a system of 
TV interacting electrons in an external potential. The idea is based on the derivation of the energy densities via 
the N — 1 -conditional probability density within the framework of the Levy-Lieb constrained search principle. 
The consequences for the development of energy functionals within the context of Density Functional Theory 
(DFT) are discussed. We propose also the possibility of going beyond the energy densities and extend this idea 
to a computational procedure where the N — 1-conditional probability is an implicit functional of the electron 
density, independently from the external potential. In principle, such a procedure paves the way for an on-the-fly 
determination of the energy functional for any system. 

PACS numbers: 02.70.Ss, 05.30.Fk, 71.15 Mb 



I. INTRODUCTION 

1. Levy-Lieb Constrained principle 

M.Levy and E.Lieb QUI] have, independently from each other, provided a general minimization principle which leads to the 
rigorous definition of the universal functional of Hohenberg and Kohn in Density Functional Theory (DFT)j3l|3l. The equation 
for the ground state energy in the Levy-Lieb (LL) formulation is: 



Egs = min 

p 



min MK + V ee \ip) + / p(r)v(r)dr 



(1) 



with Eqs the ground state energy, K the kinetic and V ee the electron-electron Coulomb operator, p(r) the one-particle electron 
density and v(r) the external potential (e.g., electron-nucleus Coulomb interaction). The meaning of EqQ]is that the minimiza- 
tion over if) is restricted to all antisymmetric wavefunctions such that p(r) = N J ip*(r, r 2 , rjv)^(r, r-2, rjv)<ir2...<ir/v, 

while the outer minimization searches over all the p's which integrate to N, number of particles. The rigorous definition of the 
universal functional of Hohenberg and Kohn follows as: 

F[p] = min {ijj\K + V ee \ip) . (2) 

Obviously, searching on the whole space of antisymmetric wavefunctions is possible only in abstract terms and becomes impos- 
sible when one tries to actually apply the LL principle and derive an explicit expression of the universal functional as a functional 
of p(r). In order to circumvent this difficulty and make it possible the derivation of a functional, one would need a formalism 
which expresses Eq|2]in terms of p, removing the explicit dependence on ip; such a formalism is reported below. 

2. The Levy-Lieb principle in terms of the (TV — 1) conditional probability density 

Let us consider the properly normalized 3./V-dimensional probability density of an TV-electron system: 

2VV>*(r,r 2 ....rjy)'0(r,r2....rjv) = 6(r, r 2 rjv) (3) 
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this can be equivalently written as: 



9(r, r 2 ....rjv) = p(r)/(r 2 , ....rjv|r) 



(4) 



where p(r) is the one particle electron density (normalized to JV) and /(r 2 , ....rjv|r) is the (JV — 1) electron conditional (w.r.t. 
r) probability density. This latter in nothing else than the probability density of finding a configuration of (JV — 1) particles after 
the position of one specific particle has been fixed 15l|6|]. In order to write the density functional in the standard notation used in 
literature, here we have identified ri with r. The LL principle of EqQ]can then be rewritten as |H1 [8], H : 



E, 



GS 



mm 

p 



( minr[p, /] 



|Vp(r) 
p(r) 



dr + p(r)v(r)dr 



with 



and 



F[ P ] 



|Vp(r)| 



p(r) 



■dr 



r[pJ] = lJ p(i 

+(JV-1) / p(r) 



|V r /(r 2 , ....rjv|r) 



R jv-i /(r 2 ,....rjv|r 

/(r 2 ,....rjy|r) 
Kv-i |r-r 2 | 



dr 2 ....dr N 
dr 2 ....dr N 



dr 



dr. 



(5) 

(6) 

(7) 
(8) 



Where R^ -1 denotes the space of configuration (r 2 , rjv). The inner minimization searches for the / which minimizes 

T[p, /]; Vp. Here we underline the fact that the above formalism does not contain approximations, i.e. the ground state identified 
in Eq. [5]is the same which solves the time-independent Schrodinger equation with the same Hamiltonian. 
The central question, is how to determine / in an efficient way and once there is a procedure for doing so, how this can be 
used in concrete terms within the DFT framework. In our previous work iflolfTTll . we adopted a physically-motivated explicit 
guess functional form for /, dependent on one free parameter, and we numerically optimized the resulting T[p, f] w.r.t. the 
single parameter. Here, we propose a radical step further, by leaving the functional form of / completely undetermined and 
(numerically) derive it within an exact quantum Monte Carlo framework. In the following part of this work we suggest two 
different but related methodologies, a) one related to the calculation of the energy density of the ground state which can then be 
used as a reference for developing analytic functionals and b) another where / can be determined as a numerical functional of p, 
independently of the external potential, and thus provide a numerically exact route to the calculation of the universal Hohenberg- 
Kohn functional. It must be taken into account that the intention of this paper is to provide a theoretical/methodological guideline 
and its practical warnings; at this stage we do not provide numerical experiments. In fact, we hope that the optimal computational 
implementation of the approach will come from a constructive discussion of the ideas reported here. 



II. ENERGY DENSITY OF THE GROUND STATE 



If we restrict ourselves to the ground state of a specific system of JV electrons with a well defined external potential, then the 
procedure of inner minimization of Eqj5](i.e. the search for / which minimizes T[p, f] , Vp) leads to f m - m = fas- 
Let us define: 



7(r) 



and 



C(r) 



|V r /(r 2 , ■■■■r iV |r)| 
R «-i /(r 2) ....rjv|r) 

/(r 2j ....r N \r] 



-dr?....dr 



N 



/r«-i |r-r 2 | 
since / m j n = fas, the explicit expression of the functional F[p] is: 



<ir 2 ....<ir 



JV 



F[P\ = / P(r) 



The term: 



1 \Vp(r)\ 2 1 

8^^ + 8 //Gs(r) + (7V ~ 1)C/Gs(r) 



l|Vp(r)| 2 1 
£(r) = 8 ~WT + 8 //Gs(r) + (N ~ 1] CfGs(r) 



dr. 



(9) 



(10) 



(11) 



(12) 



3 



is an energy density per particle expressed in terms of its kinetic f | p P (i)l — l~ ^fas( r )j an d Coulomb ((N — 1) C/ GS (r)) 
parts. 

Here, If GS (r) and C / GS (r) indicate that the quantities of Eq|9]and Eq[lO|are those calculated for the / of the ground state. 
If one knew If GS (r) and C/ GS (r) as a functional of p, this would correspond to have the universal functional of Hohenberg 
and Kohn. However, even if If GS (r) and Cf GS (r) are known only in a case-by-case situation, i.e., as functions of the position 
(and not functional of the density), the expression of Eq Q~2] would represent the energy density of the universal functional, in 
the ground state of a chosen specific system. This means that the energy density of any proposed functional in literature should 
correspond to e(r) of EqQ~2]when calculated in the ground state. At this point the key question is whether there is any rigorous 
technique which can, in practical terms (i.e., not only formally), calculate the fas and thus determine If GS (r) and C/ GS (r) @]. 
If this is the case, then for any given system - once the number of particles is fixed - one would have an explicit algorithm to 
compute the numerically exact functional of p in the ground state. In the next section, we propose that the determination of the 
minimizing fas an d of the corresponding If as (r) and C/ GS (r) can be achieved by the Ground State (GS) Path Integral (PI) 
Quantum Monte Carlo (QMC) technique. 



III. GROUND STATE PATH INTEGRAL QUANTUM MONTE CARLO FOR FERMIONS 

The Path Integral lfl2[ [T3ll ground state ifTill approach (GSPI) allows one to write the quantum partition function of a system 
of TV-particles as: 

Z = J rfR dR A fV(Ro)exp[-5(R ,Ri....RM)]V'(RAf) (13) 

here Ro = (r°, r®, r^) is a configuration of the N particles in space, equivalently Ri is another configuration and so on. In 

this way the sequence Ro Raj represents an open path of length M in the spaces of the iV-particle configurations. ip(R,o) 

and iP(Rm) is a trial wavefunction calculated at the initial and final configuration. 

Conceptually, the choice of if) is immaterial, since the evaluated quantities do not depend on it. However, technically, the choice 
of a good trial ip enhance the convergence of the method. 5(Ro, Ri ....Rm) is the action defined such that: 

cx P [-5(R ,Ri....Rm)] = <R |e- Tff |Ri)<Ri|e- TH |R 2 ) (Rm-iIc^IRm) (14) 

where r, which is formally an imaginary time, is: r = jj, with j3 formally the Boltzmann factor (the temperature has no physical 
meaning, but rather a parameter that influences the convergence efficiency of the method) and H the Hamiltonian. 
In this way, the quantum mechanical partition function is written as an integral involving a sequence of transitional probabilities 
in imaginary time r. Each of these transition probabilities can be decomposed into a kinetic part: 

(Ri|e -T ' K '|Rj+i) = ^-l^e^i^^Y (15) 

with K being the kinetic operator, and a potential part: 

<R,|e--|R !+1 ) = ^eW^,] (16) 

with V being the potential operator of the system considered. 



IV. GSPI FOR ELECTRONS AND CALCULATION OF / AND F 

In case of atoms and molecules, containing electrons (i.e. fermions) V^(R) = V ee + V ne , namely the electron-electron and the 
nucleus-electron interaction. For fermions, in the case of a real ip(R), one uses the fixed node condition: 

V fermions 

(R) = V(R) for ^(R) > (17) 
V fermlons (R) = 00 for V(R) < 0. (18) 

In case V'(R) is complex, a term is added to the free-particle part of the action. 

Since the wavefunction is defined as: Vv(R) = e^ TH ^>(Il), where yfR) is the ground state wavefunction, for r that goes to 
infinity, ip goes to the exact ground state wavefunction. Technically [12], the wavefunction is evaluated at the midpoint of the 
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path, i.e. at Rm/2- In order to proceed in the derivation of / and T in terms of the GSPI approach, we adopt the following con- 
vention: we will indicate the configuration at the midpoint of the path, Rm/2 as R*- This means that(r M / 2 ,r 2 M / 2 , vf 2 ) 

becomes (r*, r^, r^). According to EqQj] the N — 1-conditional probability density / can now be written as: 

/(i*2, ....r^|r*) = — !— J dR rfRi dRM_ 1 dRjw +1 cCRm 

V'(Ro) exp[-S(R*, R , Ri, ....Rm_ 1( Rm +1 R m )]^(Rm) (19) 

where: 

Z r , = J dR^- 1 dR dRi....dRj«_ 1 cfRM +1 ...dRMV'(Ro)exp[-S'(R»,Ro,Ri, Rm_ 1 ,Rm +1 , ....R m )]^(R m )(20) 

rfR^ -1 means that the integration is done on the whole space of configurations rj, of R* except that corresponding to 
variable r*. With this set up, / can be calculated by propagating stochastically, according to a Monte Carlo procedure, the path 
R in imaginary time r. Since the GSPI procedure, when evaluating in R„, delivers the ground state wavefunction of the system, 
the expression of / in EqfT9lcorresponds to the ground state TV — 1-conditional probability density, that is, it corresponds to / m i n 
of EqJTTJ The expression of Eq{l9]can be introduced into Eqj9]and EqJTUJ this leads to: 



and 



C(r*) = [ ^r-^W ---^. (22) 

Where now J/ min (r) = I(r*) and C/ min (r) = C(r*). Moreover, in I(r*), the gradient, V r ./(r-2, ....r^|r*) can be calculated 
analytically and thus sampled without additional computational costs in the MC sampling in configuration space; the explicit 
calculation is reported in the Appendix. The Hohenberg-Kohn functional in local form becomes: 

F[p] = \J lVp TJ/ dr* +IJ p(r*)I(r*)dr* + (N - 1) J p(r*)C(r*)dr*. (23) 

As a simple consistency check of our proposed approach, we consider the example of the homogeneous interacting electron gas. 
In this case the total Hamiltonian is K + V ee , i.e. there is no external potential v(r), thus, applying the procedure yields to the 
Local Density Approximation (LDA) approximation to the functional 0]. 



V. PRACTICAL UTILITY 

In previous work it has already appeared the idea of employing the PI approach within the framework of DFT lfl7[[T8ll . but the 
main aim there was avoiding the use of orbitals within the Kohn-Sham approach where an exchange and correlation functional, 
E X c[p], was predefined. The intention of this work, instead, is that of describing a procedure, rigorous from the conceptual and 
numerical point of view, to make it possible the numerical calculation of the exact energy density for a given system (external 
potential), and thus use this information for developing analytic functionals. The advantage of the approach used here is that 
the energy density we derive is rigorously divided in its kinetic I(r) and potential C(r) components and thus the physical 
interpretation emerges in a natural way. 

In practical terms, a possible way to use this procedure is to treat basic reference systems (e.g., single atoms, small molecules) for 
which the application of the GSPI approach is computationally feasible. This would allow for the determination of a database of 
energy densities that can be used for the development of energy functionals, and to have a novel insight into the basic physics of 
the functional in terms of each of its specific components. For instance, an accurate DFT-level description of the van der Waals 
interactions, with current functionals, for a system as simple as the helium dimer, is still an open problem |T^]. Indeed, QMC 
calculations for the helium dimer are carried on to have some understanding of such interactions with the intention of using the 
results to build better functionals on a sound physical basis (see eg. Refs. ll6lll9ll and references therein]. The approach suggested 
here would not be computationally more demanding than that of the QMC calculations of Refs. flo, [l9ll . however it would 
automatically provide the detailed (i.e. of each of the energy components) physics of the energy density of the ground state and 
thus a numerical reference in the development of energy functionals. In particular, for a given system, one may treat the problem 
for the case of interacting and for the case of not interacting electrons and calculate Ii nt (r) and C(r) for the interacting case, and 
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Inint( r ) for the non interacting case. In this way one can determine e xc (r) = Ii nt {r) + (N — l)C(r) — I n int ( r ) — / |r-r'| ^ r '' 
(where J dr' is the Hartree term), that is the exchange and correlation energy density per particle. This quantity can be 

used as a basis for developing more general expressions of e xc (r). Moreover, a newly developed functional must give the e xc (r) 
obtained by the procedure proposed here, when used for calculating the ground state of a specific system. In this context, 
relevant, long standing questions as that of the kinetic contribution to the E xc 112011 or the problem of how to extend the LL 
formulation to the kinetic energy density 112 ill could be now addressed in a more robust way. From the numerical point of view, 
the major advantage is that every time a GSPI QMC calculation is done for a given system, the quantities 7(r) and C(r) can be 
automatically determined at no additional cost. This means that one would automatically produce an increasingly larger database 
to be used for the development of functionals, and thus no more restricted to the uniform gas only. Furthermore, in this context, 
recently developed approaches within the framework of the so-called kernel-based Machine Learning (ML), propose a training 
strategy for systematically determine a numerical functional [22] (though for the moment only for a simple proof-of-concept test 
case). ML is a powerful tool for finding patterns in high-dimensional data and, applied to our strategy, would mean using GSPI 
exact results for non trivial cases to refine (train) a numerical expression of the functional. It is remarkable that the flexibility of 
ML allows for the insertion in the functional form of as much physical intuition as felt necessary (e.g., by imposing known exact 
analytic constraints), in order to reduce the dimensionality parameter space in which the numerical optimization is performed. 



VI. WARNINGS 



It must be noticed, that the use of / in QMC procedure, could lead to results characterized by a large variance. In principle, 
one may derive e(r) of Eq[T2]in more efficient manners without passing through the calculation of /, however the separation of 

the kinetic functional in the two terms: ; ^( r )^ would be no more straightforward and thus the detailed understanding of 

the physics related to the interplay between these two terms may be lost. While this is not relevant for practical applications, 
it may be relevant for the understanding of the basic physics and thus for the construction of analytic functionals. In general, 
without invoking /, in the standard GSPI procedure, C(r) can be efficiently calculated in the following way: the density is 
determined as the number of electron visiting some volume elements from the middle time slice, it follows that C (r) is the 
average electron-electron potential in that volume elements. However, it is not clear how one must then deal with the kinetic 
part (which is of major concern in this paper) although, as underlined above, it cannot be excluded that it may exists a more 
efficient way that does not directly involve /. In case such a procedure is possible, our basic idea of determining each term of 
the energy density via GSPI QMC remains valid and we would gain in computational efficiency. At the current stage this aspect 
goes beyond the aim of this paper. 



VII. UNIVERSAL FUNCTIONAL 



In the procedure discussed in the previous sections, we must restrict the calculations to cases where v(r) = V ne (r), is specified 
and then explicitly used in the QMC calculation of the transitional probabilities: 

(Ri|e-^|R i+1 ) = __i__ e -i[Vaa(R i )+Vna(R i ) + V ee( R i+1 ) + V„ e (R i+1 )]. (24) 



( 27rr )3A72 

For this reason one cannot determine F[p] as a universal functional of p independently from v(r). Here we propose to modify 
the GSPI approach so that the resulting / is a functional of p only, independently from v(r). 
To this aim we rewrite the transitional probability for the potential part as: 

1 



(Hile-^lRi+i) 



(2ttt) 3JV / 2 



r[V e e(Ri)+V ee (Ri+i)] 



(25) 



i.e., considering only the electron-electron interaction. The transitional probability for the kinetic part remains the same as above 
(Eq. H5l . Next, we can calculate /, and thus I(r) and C(r) as in Eq(l9l Eql2"Tl and Eql22l but with a sampling restricted to a 
trial ptriai- This means sampling the R^s in the configuration space with the constrain that the one particle density is ptriai, for 
more technical details about sampling at given density see note in Ref. 12311 . In this case the QMC procedure assures that the 
principle: minj T[p, /], is achieved in the sense that the resulting T is that of "ground state" at the fixed ptriai- Note that at this 
stage the external potential is not invoked and it is actually absent; in practice what we have done is to find the / of ground state 
of a gas of electron with some artificially forced electron density. From the obtained / we can now calculate the corresponding 



I(r) = 
p(r): 



Iptriai ( r ) an d C( r ) = C Ptrial (r). These quantities are taken as a first guess to write an energy functional for a generic 

'1 \Vp(ry 2 



E[p] = J P(r) 



P(r) 



-dr + -I Ptrial (r) + (N — l)C ptrial (r) + „(r) 



dr. 



(26) 
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Next, we use Eq|26] for a minimization w.r.t. p and obtain a new p = pl ut , different from ptriai because in Eqf26] the 
effect of v(r) is explicitly included during the energy minimization. At this point one can use p] mt as a new trial den- 
sity, repeat the QMC procedure, that is we search the ground state of a gas of electrons with artificially forced elec- 
tron density p\ ut , this will lead to a new / and in turn to a new J(r) and C(r) and thus we can have a new guess for 

E[p] : J p(r) | ^ V p P [r)2 dr + ^I p i t (r) + (N — l)Cy ( (r) + v(r) dr. As above we can then use the expression E[p] (again 

for a generic p(r)) for a minimization w.r.t. p(r) and obtain as a result a new p\ ut and repeat the procedure until p l out = p l J^l 
with some accuracy. Of course, convergence must be proven and intuition suggests to start from some "reasonable" ptriai- 
However, beyond the several problems that a "realistic implementation" would imply, this procedure would have at least some 
conceptual benefits; this is a real space procedure which does not require neither orbitals nor a predefinition of the functional as 
it is instead the case for the Kohn-Sham approach, but above all, the procedure has the potential to deliver the energy functional 
on-the-fiy during the calculation. Since this approach is valid for any system, independently from v(r), this is an implicit way 
to define in iterative manner the universal functional of Hohenberg and Kohn. Despite the involved computational costs are not 
clear yet, this idea may in principle offer a complementary way for using QMC in the perspective of DFT and perhaps a path to 
find a compromise between the high accuracy and computational costs of QMC and the low accuracy and computational costs of 
DFT. In any case the optimization of the computational aspects of this idea represent an interesting challenge for future research. 



VIII. CONCLUSIONS 



We have proposed a theoretical/conceptual protocol based on using the Ground State Path Integral Quantum Monte Carlo 
technique in the context of DFT for the determination of the energy functional. We have shown that the method can be certainly 
used to calculate automatically the energy density of the ground state in terms similar to those of the energy functional in 
DFT. This allows for using the results for building a database for the development and control of energy functionals beyond the 
standard case of the uniform electron gas as a reference. A second possibility to employ GSPI QMC is that of using the procedure 
as an intermediate step in an iterative loop to determine the density of ground state within the Levy-Lieb energy functional 
minimization procedure: In simple terms the functional is determined on-the-fiy during the minimization procedure. Being valid 
for any external potential, this procedure implicitly defines, in numerical terms, the universal functional of Hohenberg and Kohn. 
Despite the fact that the computational costs may turn out to be rather high, the procedure may open a way to find a compromise 
between the accuracy of QMC and the feasibility of DFT. As underlined above, next challenge would be that of searching for the 
most efficient computational implementation of the idea and compare its performance with that of standard methods. However, 
even if it will turn out to be computationally less convenient, this protocol would always assure the access to basic physical 
information that can be then employed as a complementary knowledge in the development of new, physically sound, energy 
functional; this in summary is the message of this paper. 
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Appendix 

Here we report the analytic calculations of V r » [/(rj, ....r^|r*)]. Let us define 



with 



It follows that: 



/(r*,....r^|r*) = J- x W(r*, r*, ...v* N ) (27) 



W(r*, r-2, ...r* N ) = J dR dRi c?Rm _ 1 c?Rm +1 c?Rm 

^(R )exp[-S(R*,R ,Ri,....RM !,Rm +1 R m )]^(Rm) (28) 



V r 4/(r*,....r^|r*)] = - *_W( r *, r* 2 , ...r* N ) x V r ,Z r , + -±- x V r .W(r* ,r* 2 , ...r* N ) (29) 
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with 



V r « Z r * — J dR^ ^cIRqcIRi cZR,M._^dR,_M c?Rm 

^(Ro)V r . (exp[-5(R», R , Ri, ...JRu._ 1 ,B.^. +1 R M )]) ^(Rm)- (30) 

In, V r * ^exp[— 5(R», Ro, Ri, ....Rm_[,Rm + ] R^-)]V the terms which are involved into the derivation are only those 

with R*. This means that the interesting quantity to calculate is: V r * (R.M_i|e~ Tff |R* \/ H*\e~ TH |Rm +1 



terms of exp[— S 1 ] factorize. This gives: 



the other 



Rm._iI<5 rff |R*WR„|e Ttl |Rm_|_ 1 



-tHi 



( 27rr )3AT/2 



R M -R* \ 2 /R*-R 



For the kinetic part one has: 



that is: 



V,. 



( 27rr )3iV/2 



R. — R. # \ ^ X R- * — R- 



-f[V(R„_ 1 )+V(R»)] -5 [V(R,)+V(Rm +1 )] 

e 2 e 2 +± 



R. — R^ \ ^ / * — 



(31) 



(27rr) 3JV / 2 



r * M 1 . . , . Mil.- 

r [( r --! _ r *)_ ( r * _ r -+!)] (32) 



( 27rr )3JV/2 



R- — R- * \ ^ / R-* — R. 



r(r 2 +1 — r 2 1 )x 



( 27rr )3iV/2 



R> — R- * \ ^ / R- * — R- 



(33) 



For the potential part: 



( 27rr )3A72 



V(R.* )+V^(R. m. _ 1 ) J -f y(R,) + V(R M 



-TV r .y(R») x 



I ^-5 ^(r,)+\/(Rm +1 )J ^-f ( v(R»)+y(R # _ i : 



Mil M T , , 

( r -+! _ r --i) _ rV r .F(R» 



( 27rT )3W/2 

It follows: 

V r * Z v * — j ' dR^ ^c^RgfiRi (iR_M _^(iR_M c?R^f 

x^(Ro) exp[-S(R«, R , Ri, ....Rm _ 1; Rj« +1 R Af )]i/'(Rm) 

It follows also that: 

V T *W{r*,r* 2 ,...r* N ) = J dR dR 1 ....dR,M_ 1 dKM +1 dR M (r^ +1 - r^- 1 ) - r V r .y(R« 

xi(j(R ) exp[— 5(R*, R , Ri, ....Rm _ v Rm +1 R A/ )]^(R M ) 



(34) 



(35) 



(36) 



The term (r^~ +1 — r"^ -1 ) can be calculated without any additional computational cost during the sampling, however this 
is true also for the term V r »V r (R„). In fact the explicit form of 1/(R„) is known and thus its gradient can be calculated 
analytically and sampled without additional computational costs. If one uses the expression of Eql35l and Eq|36]into I(r*) 
(C(r*) is straightforward) obtains the exact form of the Hohenberg-Kohn functional. 
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